Do invasive herbivore species threaten ongoing Aboriginal cultural connection to Country?

Acknowledgments

We thank Chris Griffiths, Pam O’Brien, Tristan Ricketson, Jess Ward-Jones, Ian Pulsford, David Freudenberger, David Berman, and Joanne Lenehan.

This document was based on the Draft Ecology and Conservation Modelling Preregistration Template by Chris Jones and Elliot Gould. We omitted a number of their proposed sections for which we had no answer at this stage - look out for their full template and follow the project at xxxx.

1. Problem Statement

1.1 Study context and Purpose

1.1.1. Key Stakeholders

This research project pre-registration proposes an investigation on the Country of the Monaro Ngarigo people, Traditional Custodians of land including the southern section of Kosciuszko National Park, NSW, Australia. The application of the concept of significant cultural values described here is informed by people who speak for that Country and their trusted contacts.
The ongoing custodianship of the Monaro Ngarigo people is acknowledged by NSW Department of Planning, Industry and Environment’s National Parks and Wildlife Service (NPWS) through a Memorandum of Understanding for joint management of the southern section of Kosciuszko National Park [1]. An earlier version of this proposal was presented to the Executive Committee representing the Aboriginal Community of the Southen Snowy Mountains (the group as named in the MoU) for consultation, and they should be encouraged to remain involved in the project should it proceed.

1.1.2 The problem, and purpose of the model

Feral herbivores pose a serious threat to Australian ecosystems and native biodiversity [2], including our species and ecological communities threatened with extinction or collapse. There is far less attention to how those herbivores may threaten significant cultural values of landscape for Aboriginal Australians.

The landscape and riverscape of the lower Snowy River Valley has by many published and anecdotal accounts been increasingly impacted by feral herbivores [3]. Recent dung counts from a set of exclosures established in the 1980s to distinguish between the impacts of rabbits and kangaroos instead tell the story of the increasing dominance of horse and deer [4]. Even after a period of drought followed by fierce wildfire, feral horse numbers in Kosciuszko National Park were estimated at 14,000 in Spring, 2020 [5]. Deer populations have also been expanding through the area and Herbivore signs are everywhere from the valley floor to the stony ridges.

This study will examine the relationship between vertebrate herbivores - feral and native - and the regeneration of kurrajong (Brachychiton populneus), a tree species that exemplifies the living connection of Monaro Ngarigo people to Country. A focus on one particular cultural value risks missing the point of the deep and encompassing spiritual connection of those people to Country as a living whole [6], but at least this particular example of a tree species may encourage readers to think of living cultural connection to landscape rather than a historical presence evoked by a focus on physical artifacts.

The kurrajong is what is referred to as a ‘culturally significant species’ (sensu [7]), an important resource providing food (seeds, new leaf buds), fibre (fish nets, string), gum, and the possibility of water in hard times whose value for Aboriginal Australians has been documented for a long time (e.g., [8]). Kurrajong is widely distributed in south-eastern Australia and is not considered at strong risk of extinction, but global scarcity is not a criterion for cultural significance [7]. For Monaro Ngarigo people kurrajong probably signified no less than food and fibre, and probably more. The canopy of large individuals provide remarkable cooler and moister conditions than the surrounding rain-shadow woodlands (authors observation), and kurrajong are reported to have been placed and encouraged in certain locations, either through translocation of seedlings, or dispersal of seed [9]. That practice may explain in part the current day distribution of large old kurrajong dotted through the lower Snowy River valley including up into the drier slopes.

In order to manage the herbivore threat to a particular value, one would ideally know the relationship between herbivore population density and the impact on that value [11]. These density-impact functions could take a range of forms (e.g., [11]] and knowing which is applicable has enormous practical implications for managers. A highly sensitive function might suggest that only reducing herbivore populations to very low density will adequately protect the value (e.g., [12] for rabbits), whereas for a relatively insensitive function, avoiding irruptive population growth is adequate. An understanding of the density-impact function therefore allows a decision-maker to design an intervention likely to achieve the required benefit with the minimum financial cost [10] and minimum destruction of herbivores. This latter point is vitally important as the removal or destruction of feral horses is controversial.

A density-impact function requires that the value that is being impacted be clearly defined. The problem that high densities of feral herbivores may pose for kurrajong is persistent failure to establish a healthy new cohort of saplings. Mature old kurrajong are venerated, but kurrajong offer useful resources from sapling age: supple fibre could be easier to extract from the bark of saplings than from older individuals, and fresh leaf buds could be easily plucked for eating from accessible low foliage. Thus, the demographic structure of kurrajong at a point in time, and evidence for the survival and growth of emerging generations (the cultural connection through the value being maintained), represent the scope of an operational definition of the value for this study.

This observational study will examine the relationship between vertebrate herbivore activity and the demographic structure and seedling survival of kurrajong in the lower Snowy River Valley. The study will be limited to locations where veteran kurrajong are found on the basis that those locations represent the most likely site where seedling establishment and progression through sapling phase to maturity would have occurred. If a strong association between herbivore activity and demographic or seedling frequency and condition is found, we will explore what levels of herbivore activity (a proxy for density) present a threat to ongoing population stability.

1.1.3 Research questions and hypotheses

Is kurrajong demographic structure related to feral herbivore activity? If feral herbivores affect demographic structure, then all else being equal, higher herbivore activity will reliably ‘predict’ sites where veteran kurrajong are accompanied by few seedlings and saplings, and those saplings will show signs of having been repeatedly browsed or damaged.

Kurrajong recruits readily from seed and is robust to transplanting at young seedling age. In some circumstances it has proved capable of prolific if not invasive recruitment under established canopy of other species e.g., [13]. Therefore seedling recruitment appears likely to be a dynamic response variable in this context.

Veteran kurrajong are to be found on the valley floor and hillslopes in the lower Snowy Valley corridor up to around 750 m altitude (Fig. 1). Considerable variation in soil depth, nutrition, and water availability is expected across that natural gradient from mild slope –> valley floor –> exposed harsh slope, and that variation would be expected to influence inherent germination and survival probabilities for kurrajong seedlings in a hypothetical herbivore free setting, and also influence herbivore activity (better resource availability promoting foraging opportunity, therefore higher activity).

Figure 2. Basic sampling concept, based around veteran kurrajong on hillslopes and the valley floor.

A main effect of sample location on kurrajong demographic structure and seedling survival is anticipated, and a main effect of herbivore activity. An interaction is expected whereby relatively less herbivore activity in lower resource locations has a disproportionately high impact on kurrajong demographic structure and seedling survival.

1.1.3 Development of response variables

One measure of interest is population demographic structure, the distribution of age classes at sampling locations, which may relate to contemporary and historic patterns of herbivore activity.

However, as each sampling unit is intended to be based on a veteran individual, and at the intended sample size, we are unlikely to observe sufficient individuals to generate robust statistics of distribution per plot.

Response variables for modelling will likely be fashioned around the density of juvenile individuals around the veteran individual. If we imagine the following classes — new seedling (first season) –> established seedling (1-2 years) –> sapling (2-5 years) –> pre-reproductive (5-\(x\) years) –> mature reproductive (\(x\)–\(y\) years) –> veteran (\(>y\) years) — we might be particularly interested in the first three categories.

Seedling density

The number of “seedlings” in a quadrat is one simple candidate response variable. However, very small plants could be newly germinated (e.g., Fig. 3 - recently germinated individual with cotyledons) or an older ‘seedling’ indicating various seasons of growth and repeated browsing (Fig. 4).

Seedling / young sapling intactness

The density or number of seedlings could be augmented with evidence about past browsing damage to stems. Stem diameters above and below ground is one way of looking for evidence of severe past browsing, in addition to the more obvious above ground signs such as in Fig. 4.

Condition and age of older saplings and mature individuals

Size, and browsing and rubbing damage indices could be collected for older saplings and trees. Sampling tree cores for dendrochronological analysis would be ideal to be able to calibrate tree size with approximate age. If ageing can be done convincingly, it might be possible to explore whether the apparent cohorts gap between older saplings and mature trees corresponds to environmental conditions or estimates of historical grazing pressure.

Ground layer structure and cover

We could try to include some assessment of the cover or biomass of an understorey indicator species attractive to herbivores, like the native pea Swainsona, the vanilla lily Arthropodium or kangaroo grass (Themeda).

2.0 Conceptual Model and Study Design

2.1 Conceptual Model Representation and Elicitation Method

Our conceptualisation of the problem is informed by a field trip undertaken by the authors to the southern part of Koscsiuszko National Park in December, 2020. In addition, author RS has a long association with the landscape recreationally, working as a professional river guide, and working as a contractor to the parks management agency. The proposal, laying out the model has been referred to NPWS (the park manager) and to the representatives of the Monaro-Ngarigo clan on the Executive Advisory Committee.

2.2.1 Explain Critical Conceptual Design Decisions

One of the unavoidable complications of studying contemporary vertebrate herbivore pressure is the legacy of damage from the livestock grazing era (Pulsford, pers. comm). The current day impacts of feral herbivores is superimposed on whatever legacies of the era when sheep and cattle where driven into the system. Whatever legacy may remain from that period is unavoidable, but will be important to bear in mind in project design and analysis. Our intention is to try replicate the design further upstream beyond around Black Jack Gully (see Fig. 5) where apparently the evidence of stock droving is far less.

Figure 5. Case study area, the Snowy River Valley in the southern Kosciuszko National Park. The image below indicates three sections of the Snowy corridor that might have different patterns of use. We propose to sample in each section in the hope of teasing apart legacy effects of livestock droving from contemporary effects of feral herbivores.

We are principally interested in herbivore density but in practice will measure dung deposits as an index of herbivore activity. Dung deposits across a range of vertebrate herbivore species in the Australian environment are typically assumed to last up to 12 months, thus they may notionally integrate about that period of time. We will rest on the most appropriate estimates of dung residence available for each species rather than seeking to establish the decomposition profiles afresh.

2.2.2 Model assumptions and uncertainties

Persistent soil moisture deficit or drought is one important exogenous variable that is beyond our control in this study. Soil moisture availability has a strong influence on ground layer vegetation, seedling survival, herbivore foraging patterns, herbivore population density, and herbivore density impact functions. For example, during the drought leading up to the end of 2019 – an important antecedent condition contributing to the mega-fires of 2019/20 – the impact of increasingly desperate feral horse and deer populations was at unprecedented levels [3].

2.3 Predictor Variables

Herbivore density

The density of feral vertebrate herbivores is the main predictor variable of interest, but in fact the dung count data we will gather is an index of activity, from which animal density estimates may be calculable. We will count dung deposits for horse, deer, macropods, rabbits and any other identifiable type.

Whilst a long-run ecological study of herbivore exclosure on native vegetation within our study landscape uses dung dry mass as an index[4], we expect to use dung counts, such as described by [14] and [15]. The advantage of counting is that we can avoid needing to transport large quantities of dung between remote sites and then back to a vehicle or canoe.

Elevation

Kurrajong are thought to ‘top out’ at around 750 m ASL, a limit which is encountered within the lower Snowy River valley, so the elevation of sample locations will be recorded.

Soil characteristics

Soil variables influence the availability of resources that support plant growth. We will need to understand the nature of the difference in the different valley sampling positions but we do not yet have a definitive plan for how to do this, and what mix of field and modelled variables might be usefully combined.

2.5 Sampling Plan

Veteran kurrajong are to be found on the valley floor and hillslopes in the lower Snowy Valley corridor up to around 750 m altitude. Considerable variation in soil depth, nutrition, and water availability is expected across that natural gradient from mild slope –> valley floor –> exposed harsh slope, and we will sample in each because we are interested in whether the activity-impact on demographic structure and seedling survival relationship is consistent or not (see Fig. 1).

This study aims to explore whether the contemporary pattern of herbivore activity may threaten kurrajong regeneration and future population replacement, but there are other important factors in our study area that may exert an influence.

The study area has been exposed to various kinds of vertebrate herbivore pressure over more than 100 years, including management of cattle and sheep. The legacy effects of historical phases of managed and wild herbivore access are unavoidable, but we may shed some light on this by replicating the sample design in various positions through the lower Snowy River corridor (see Fig. 5). Certain indices of intactness suggest that early droving activity was concentrated below Black Jack’s Creek, and that below where the river meets the road additional impacts from car-based access are possible.

2.6 Data collection procedures

Data collection would occur ideally in a Spring season for ease of identification of vascular plants co-occurring with kurrajong.

Data collection will be based on a 50 m x 10 m quadrat centred on a veteran kurrajong (see Fig. 6). The quadrat size was largely determined with dung detection in mind, following conversation with Dave Berman, and the work of Joanne Lenehan [15] at Guy Fawkes National Park.

Figure 6. Anticipated quadrat design for assessment of kurrajong recruitment, vegetation structure and herbivore activity covariates. Circles represent kurrajong of different ages, and blobs are dung.

Identification of veteran kurrajong (sampling locations)
We will select veteran kurrajong individuals in as objective a manner as we can, though the practical difficulties in reaching some parts of the study area will impose some limits.

Initially, we will try to identify prospective locations based on aerial imagery, at least for the hillslope sites, as kurrajong canopy tends to stand out bright green from the subdued canopy tones of the dominant species yellow box, white box and native pine. That process will be augmented to achieve the desired sample size by stopping the vehicle at predetermined locations and sampling the nearest veteran individual in the desired hillslope / valley floor location.

Rationale for excluding prospective quadrat locations
Outside of any restrictions imposed by the NPWS who would issue the research permit, or field safety risk mitigation concerns, we don’t foresee a reason to exclude pre-selected locations or individuals, unless there turns out to be no veteran kurrajong present.

quadrat position
The quadrat would - as a starting proposition - be oriented to align with the principal direction of resource movement, i.e., aligned to the slope on the hill sites, and to the stream on the valley floor. The idea being that seeds may disproportionately accumulate downstream or down-slope. However, given that the hillslopes are quite steep, if deer and horse are thought to approach likely foraging opportunities along the contour then there may be an argument to align the hillslope quadrats perpendicular to the slope.

2.4.2 Data Processing and Preparation

Numeric predictors will be centred and scaled to facilitate comparison of coefficients (effect sizes).

Assuming soil and vegetation structural attributes were measured, these may be reduced to one or two summary variables via Principal Component Analysis to test their utility in the model compared with individual attributes.

The seedling damage variable will most likely require coding as an ordinal variable where 0 = no damage… n = seedling presumed dead. These categorisations need to be developed in advance of data collection, or at least the set of instructions for how many classes and on what basis they will be delineated, to avoid the perception of convenience definitions in service of a particular result. This activity requires a quick literature review.

2.4.3 Data Exploration or preliminary data analyses

The age-size distribution of kurrajong at each location will be at least visually explored, though it seems likely that the modelled response will only be possible with seedling number or number x condition responses.

The continuous covariate data will be visually explored as a function of valley position and river section as part of familiarisation and characterisation of the sample locations.

2.4.4 Data evaluation, exclusion, and missing data

The form and distribution of response and predictor variables will be visualised to help confirm their reliability for the analytical purpose, and compared to available reference values from other studies (e.g., dung data with Lenehan [[15]), floristic data with Ward-Jones et al [4].

Any perceived issues of data reliability or credibility will be documented in the project codebook.

Where data or sample locations need to be excluded, the effect of their exclusion on the interpretation will be noted in results, with the alternative analysis available as supplementary information.

3.0 Analysis Plan

3.1 Model Class, Framework and Approach

3.1.1 Operationalising Model Variables

We will model the mean number of seedlings / seedlings + saplings encountered per unit area, assumed to be drawn from a Poisson distribution.

Seedling damage will be recorded as an ordinal variable according to a predefined set of attributes. The modeled quantity is the log-odds of being in an equal or lesser damage class.

Analysis of the demographic structure may be exploratory only as the number of individuals per quadrat may be insufficient to characterise the distribution function of numbers against size/age.

3.1.2 Model Structure

functional forms of interactions

The main interaction we foresee is the possible interaction of herbivore activity measure and river section / sampling location. In principle

discretisation of continuous variables

The main discretisation that we foresee is the recording of browser damage to seedlings (notionally a continuum ) into discrete damage classes.

3.1.3 Model Class / Family

Specify which family of statistical distributions you will use in your model, and describe any transformations, or link functions.

• Include in your rationale for selection, detail about which variables the model outputs are sensitive to, what aspects of their behaviour are important, and any associated spatial or temporal dimensions in sampling.

3.2 Model Structure and Parameter Estimation

3.2.1 Parameter/Structure estimation technique

Candidate model structures will be guided largely by the ecology of the system. We will report on the models described in this document even where candidate predictor variables do not feature in the most parsimonious statistical model (though that too will be reported).

This being our first study pre-registration, it is likely that we do not yet have well honed powers of foresight, therefore it is also likely that the act of analysis will reveal some erroneous assumptions outlined at this juncture.

3.2.2 Estimation performance criteria

Important model features or coefficients will be identified by the magnitude and precision of their effect on the predicted number of kurrajong seedlings in good physical condition.

3.3 Model assumptions and uncertainties

• Specify all assumptions and key uncertainties in the formal model.

• Describe what gaps exist between the model conception, and the real-world problem, what biases might this introduce and how might this impact any interpretation of the model outputs, and

3.4 Specify formal model(s)

The seedling encounter rate per quadrat \(R\) will be assumed to be drawn from a Poisson distribution with a grand mean \(\lambda\).

\[ R_{seedlings} \sim dpois(\lambda)\\ log(\lambda) = soil + AI_{deer} + AI_{horse} + I_{macropod} + S * P + \epsilon_i\]

where \(AI\) is the activity index for the various herbivore species, \(S\) is the river section and \(P\) is the valley sampling position.

Seedling damage class \(D\), an ordinal variable with damage classes \(j\). The odds of a seedling being in a damage class less than or equal to than any given class is by the following equation:

\[ \frac{P(D\leq j)}{P(D>j)} \]

The odds are modelled using the log-odds or logit function:

\[ log\frac{P(D\leq j)}{P(D>j)} = logit(P(D\leq j)) \]

The unit of observation for this model is seedling \(i\), which are nested within sample quadrats \(q\). This corresponds to a generalized linear mixed model structure.

\[ logit(P(D\leq j)) = AI_{horse} + AI_{deer} + AI_{macropod} + S * P + Elev + \epsilon_q \]

4.0 Model Calibration, Fitting & Checking

4.1 Model calibration and validation scheme

Internal validation

The model of seedling density will be validated with leave-one-out cross validation, where the model’s skill in predicting the seedling density in each held-out quadrat is summarised for the sample \(n-1\) iterations.

The model validation for the seedling condition model remains to be determined.

External validation

Depending on the time and resource availability we could select an additional set of new sites and test the predictive performance on the seedling densities and seedling condition observed at those sites.

The ability to manipulate herbivore densities through structures like exclosures might allow an additional validation opportunity in the future, but that is beyond the scope of this project.

4.1.2 Describe calibration/validation data

At this stage we have not considered how this might be achieved for our study.

4.2 Model Checking

Quantitative checking and evaluation are not simple for this class of model. The goodness of fit tests applied to this mixed model of seedling condition will be checked following [16] unless some extension can be identified using simulated quantile error residuals via the DHARMa package [17].

4.3.1 Quantitative model checking

Observed data will be plotted over fitted / adjusted model estimates to identify characterise the nature of any discrepancy and or failure to represent our study system.

4.3.2 Qualitative model checking

The model predictions will be presented as accessible visualisations with easy-to-interpret scales to the community of scientists and managers.

4.3.3 Assumption violation checks

At this stage we have not considered how violation checks might be done.

5.0 Further information

An accessible summary of this case study is available at nespthreatenededspecies.edu.au [XXXXX]

References

[1]

“Southern snowy mountains aboriginal community MoU. Kosciuszko national park plan of management,” 2016. https://www.environment.nsw.gov.au/-/media/OEH/Corporate-Site/Documents/Parks-reserves-and-protected-areas/Aboriginal-joint-management/southern-snowy-mountains-aboriginal-community-memorandum-of-understanding.pdf (accessed Feb. 17, 2021).

[2]

S. G. Kearney et al., “The threats to australia’s imperilled species and implications for a national conservation response,” Pacific Conservation Biology, vol. 25, no. 3, pp. 231–244, 2019, doi: 10.1071/PC18024.

[3]

J. Blay, “The curse of the wild horses: De-romanticising a plague of feral wildlife (to save australia’s unique kosciuszko national park and its wildlife).”

[4]

J. Ward-Jones, I. Pulsford, R. Thackway, D. Bishwokarma, and D. Freudenberger, “Impacts of feral horses and deer on an endangered woodland of kosciuszko national park,” Ecological Management & Restoration, vol. 20, no. 1, pp. 37–46, 2019, doi: 10.1111/emr.12353.

[5]

U. Malone, “New survey supports need to cull more of kosciuszko’s wild horses, state government says,” Jan. 2021, Accessed: May 03, 2021. [Online]. Available: https://www.abc.net.au/news/2021-01-13/wild-horse-population-in-kosciuszko-slashed-by-a-quarter/13053240.

[6]

Kosciuszko Aboriginal Working Group, “Yerribie / dhirrayn,” in Kosciuszko national park plan of management: 2006, Sydney South, N.S.W.: Dept. of Environment; Conservation NSW, 2006.

[7]

I. R. G. of the National Environmental Science Program’s Threatened Species Recovery Hub, “A case for culturally significant species: Submission to the 2020 independent review of the environment protection and biodiversity conservation act 1999.” 2020.

[8]

J. H. Maiden, The useful native plants of australia:(including tasmania). Turner; Henderson, 1889.

[9]

J. L. Silcock, “Aboriginal translocations: The intentional propagation and dispersal of plants in aboriginal australia,” etbi, vol. 38, no. 3, pp. 390–405, Sep. 2018, doi: 10.2993/0278-0771-38.3.390.

[10]

H. Yokomizo, H. P. Possingham, M. B. Thomas, and Y. M. Buckley, “Managing the impact of invasive species: The value of knowing the density–impact curve,” Ecological Applications, vol. 19, no. 2, pp. 376–386, 2009, doi: 10.1890/08-0442.1.

[11]

G. L. Norbury, R. P. Pech, A. E. Byrom, and J. Innes, “Density-impact functions for terrestrial vertebrate pests and indigenous biota: Guidelines for conservation managers,” Biological Conservation, vol. 191, pp. 409–420, Nov. 2015, doi: 10.1016/j.biocon.2015.07.031.

[12]

G. Mutze, B. Cooke, and S. Jennings, “Estimating density-dependent impacts of european rabbits on australian tree and shrub populations,” Australian Journal of Botany, vol. 64, no. 2, p. 142, 2016, doi: 10.1071/BT15208.

[13]

M. Buist, C. J. Yates, and P. G. Ladd, “Ecological characteristics of brachychiton populneus (sterculiaceae) (kurrajong) in relation to the invasion of urban bushland in south-western australia,” Austral Ecology, vol. 25, no. 5, pp. 487–496, 2000.

[14]

B. Walters, “Report on the feasibility of assessing feral horse densities in the Australian Alps Parks using a strip-transect method,” Hurstbridge, Vic, 1996.

[15]

J. R. Lenehan, “Ecological impacts of feral horses in grassy woodland and open-forest gorge country in a temperate-subtropical wilderness,” PhD thesis, University of New England, 2011.

[16]

D. W. Hosmer, S. Lemeshow, and R. X. Sturdivant, Applied Logistic Regression, Third. John Wiley & Sons, Ltd, 2013.

[17]

F. Hartig, DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. 2017.